###TWAS结果展示,属于疾病类型
library(ggrepel)
library(ggplot2)
library(dplyr)						#第一部分的代码用不了了，因为这个结果是meta分析的，meta只能用于不同研究的荟萃分析。我记得case是一样的。
; setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201204")
; ASD=read.csv("meta.p.TWAS.ASD.statis.csv",header=T)	#这是阿尔兹海默症的数据
; BD=read.csv("meta.p.TWAS.BD.statis.csv",header=T)
; Depr=read.csv("meta.p.TWAS.Depr.statis.csv",header=T)
; SCZ.BD=read.csv("meta.p.TWAS.SCZ.BD.statis.csv",header=T)
; Depr=Depr[,-6]
; ASD$group.source="Alzheimer"
; BD$group.source="BD"
; Depr$group.source="Depr"
; SCZ.BD$group.source="SCZ.BD"
; rt=rbind(ASD[1,],BD[1,],Depr[1,],SCZ.BD[1,])
; rt$logpvalue=-log10(rt$P.value)
; rt$counts=c(296,645,454,2417)
; rt$group=ifelse(rt$P.value<0.05,"sig","nosig")
; rt=rt[order(rt$P.value),]
; rt$num=4:1
; rt$logcounts=log2(rt$counts)

; ggplot(data=rt,aes(y=OR,x=num))+geom_errorbar(aes(ymin=lower,ymax=upper,color=P.value),width=0.1)+geom_point(aes(size=0.5,color=P.value))+geom_text_repel(aes(label = counts))+
    ; scale_color_gradient(low = "#E65B35", high = "#3C5488")+coord_flip()+geom_hline(yintercept = 1,color="red",linetype="dashed")+scale_x_continuous(breaks = rt$num,labels = rt$group.source)+
    ; theme_bw()+ylab("OR")+scale_y_continuous(breaks=c(seq(0,2,0.5)),limits = c(0,2))+
    ; theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank(),axis.text.x = element_text(size = 15,color="black"),axis.text.y = element_text(size = 15,color="black"))

						#这是第二部分，只取random1做部分结果的展示

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201204")
ASD=read.csv("TWAS.ASD.random.con.1.statis.csv",header=T)[6,]	#这是阿尔兹海默症的数据,缩写应该为AD
BD=read.csv("TWAS.BD.random.con.1.statis.csv",header=T)[6,]
Depr=read.csv("TWAS.Depr.random.con.1.statis.csv",header=T)[6,]
SCZ.BD=read.csv("TWAS.SCZ.BD.random.con.1.statis.csv",header=T)[6,]
ASD$group.source="Alzheimer"
BD$group.source="BD"
Depr$group.source="Depr"
SCZ.BD$group.source="SCZ.BD"
rt=rbind(ASD,BD,Depr,SCZ.BD)
rt$logpvalue=-log10(rt$P.value)
rt=rt[order(rt$P.value),]
rt$num=4:1
																#散点图好看些
library(ggplot2)
library(ggrepel)
ggplot(rt,aes(logpvalue,OR))+geom_point(aes(size=3,color=logpvalue))+
  scale_color_gradient(low = "#3C5488", high = "#E64B35")+
  theme_light(base_size = 15)+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+
  geom_text_repel(data = rt,aes(label = paste(group.source,overlap.num,sep="\n")),size = 5)+
  theme(panel.grid.minor = element_blank())+scale_x_continuous(breaks=c(seq(0,3.5,0.5)),limits=c(0,3.5))+
  xlab("log10(P-value)")

  rt$group=ifelse(rt$FDR<0.1,"significant","nosignificant")
rt$group=factor(rt$group,levels = c("significant","nosignificant"))
ggplot(rt,aes(x=num,y=overlap.num,group=group))+geom_bar(stat ="identity",width = 0.5,aes(fill=group))+
    scale_x_continuous(breaks = rt$num,labels = rt$group.source)+scale_fill_manual(breaks = c("significant","nosignificant"),values = c("#E64B35","#3C5488"))+
    geom_text(aes(label = overlap.num),position=position_dodge(width = 1),size = 5,vjust = -0.25)+coord_flip()+theme_classic(base_size = 15)
yintercept
xintercept
###染色质状态结果展示，属于调控类型
library(dplyr)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201222")
rt=read.csv("HMM.enrich.117K.rm.eQTL.support.psyASH.meta.csv",header=T)
rt1=rt[rt$NUM.state<10,]
rt2=rt[rt$NUM.state>=10,]
rt1=arrange(rt1,group,P.value)
rt2=arrange(rt2,group,P.value)
rt=rbind(rt2,rt1)
rt$num=1:15
rt$logOR=log2(rt$OR)
rt$logupper=log2(rt$upper)
rt$loglower=log2(rt$lower)
rt$group=ifelse(rt$P.value<0.05,"sig","nosig")

ggplot(data=rt,aes(y=logOR,x=num))+geom_errorbar(aes(ymin=loglower,ymax=logupper,group=group,color=group),width=0.5)+geom_point(aes(size=0.5,color=group))+
    scale_color_manual(values = alpha(c("#3C5488","#E64B35"),0.9))+coord_flip()+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_x_continuous(breaks = rt$num,labels = rt$DESCRIPTION)+
    theme_bw()+ylab("log2(OR)")+
    theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank(),axis.text.x = element_text(size = 15,color="black"),axis.text.y = element_text(size = 15,color="black"))

  
###H3k的组蛋白甲基化修饰出现的可能性与羟甲基化方向，属于调控类型
library(ggplot2)
library(dplyr)
library(tidyr)
#	807psy的结果,打算在文章中也用807的结果，因为它结果好一点，H3K4me3属于转录激活的，这个修饰状态下same的数量为438，
#	H3K3me3是抑制的，此状态的opposite为451个，same为356个

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201227.H3k.analysis/807.analysis/")
file=read.csv("807.psy.vs.healthy.H3k.statis1.csv",header=T)
file$group=ifelse(file$P.value<0.05,"significant","nonsignificant")
file=arrange(file,group,OR)
#file=file[order(file$P.value),]
file$num=1:7

ggplot(data=file,aes(y=OR,x=num))+geom_errorbar(aes(ymin=lower,ymax=upper,group=group,color=group),width=0.5)+geom_point(aes(size=0.5,color=group))+
  scale_color_manual(values = alpha(c("#3C5488","#E64B35"),0.9))+coord_flip()+geom_hline(yintercept = 1,color="red",linetype="dashed")+scale_x_continuous(breaks = file$num,labels = file$H3k.group)+
  theme_bw()+ylab("OR")+scale_y_continuous(breaks=c(seq(0.8,1.5,0.1)),limits=c(0.8,1.5))+
  theme(panel.grid.minor = element_blank(),panel.grid.major = element_blank(),axis.text.x = element_text(size = 15,color="black"),axis.text.y = element_text(size = 15,color="black"))

#柱状图
file=read.csv("807.ASH.analysis.csv",header =T)
file$group=ifelse(file$P.value<0.05,"significant","nonsignificant")
file=file[order(file$P.value,decreasing = T),]
file1=data.frame(file[6:7,1:3])
file1$num=1:2
file1=gather(file1,key = direction.group,value = counts,same.mean,opposite.mean)

ggplot(file1,aes(x=num,y=counts,,group=direction.group))+geom_bar(stat ="identity",width = 0.5,position = "dodge",aes(fill=direction.group))+
  scale_fill_manual(values = alpha(c("#3C5488","#E64B35"),1))+scale_x_continuous(breaks = file1$num,labels = file1$H3k.group)+
  geom_text(aes(label = counts),position=position_dodge(width = 1),size = 5,vjust = -0.25)+coord_flip()+theme_classic(base_size = 15)


#8544 ASH的结果
library(ggrepel)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201227.H3k.analysis/8544.analysis/")
file=read.csv("8544.ASH.analysis.csv",header = T)
file$group=ifelse(file$P.value<0.05,"significant","nonsignificant")
file=file[order(file$P.value),]
file$num=7:1
file$logpvalue=-log10(file$P.value)
#散点图，p值和same.ratio为坐标轴
ggplot(file,aes(y=same.ratio,x=logpvalue,color=group))+geom_point(size=3)+scale_color_manual(values = alpha(c("#3C5488","#E64B35"),0.9))+
  labs(y="same_ratio",x="-log10 (Pvalue)",title="H3k")+
  theme_classic(base_size = 15)+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+
  geom_text_repel(data=subset(file,file$P.value<0.05),aes(label=H3k.group),size=3, fontface="bold",force = T,box.padding = unit(0.5, "lines"),point.padding = unit(0.8, "lines"), segment.color = "black", show.legend = FALSE)

#柱状图
library(tidyr)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201227.H3k.analysis/8544.analysis/")
file=read.csv("8544.ASH.analysis.csv",header = T)
file$group=ifelse(file$P.value<0.05,"significant","nonsignificant")
file=file[order(file$P.value,decreasing = T),]
file$num=seq(1,13,2)
file1=data.frame(file[4:7,1:3],file[4:7,5:7])
file1=gather(file1,key = direction.group,value = counts,same.mean,opposite.mean)

ggplot(file1,aes(x=num,y=counts,,group=direction.group))+geom_bar(stat ="identity",width = 1,position = "dodge",aes(fill=direction.group))+
  scale_color_manual(values = alpha(c("#3C5488","#E64B35"),1))+scale_x_continuous(breaks = file1$num,labels = file1$H3k.group)+
  geom_text(aes(label = counts),position=position_dodge(width = 1),size = 5,vjust = -0.25)+coord_flip()+theme_classic(base_size = 15)


### gene set的富集，2021.01.01
library(ggplot2)
library(ggrepel)
result=read.csv("./20201231.GeneSet/GeneSet_enrichment.807psyASH.csv",header = T)
rt=result[result$con.file=="con_genotype.random1.anno",]	#只取random1做部分结果展示
rt$FDR=p.adjust(rt$p.value,method = "BH")
rt$logFDR=-log10(rt$FDR)
ggplot(rt,aes(logFDR,OR))+geom_point(aes(size=3,color=logFDR))+
  scale_color_gradient(low = "#3C5488", high = "#E64B35")+
  theme_light(base_size = 15)+geom_vline(xintercept = -log10(0.1),color="red",linetype="dashed")+
  geom_text_repel(data = subset(rt,FDR<0.1),aes(label = paste(geneset,case_overlap,sep="\n")),size = 5)+
  theme(panel.grid.minor = element_blank())+scale_x_continuous(breaks=c(seq(0,3.5,0.5)),limits=c(0,3.5))+
  scale_y_continuous(breaks=c(seq(0,3,0.25)),limits=c(0,3))+xlab("log10(FDR)")
  
###TF analysis 采用807预测的TF，然后分析其与羟甲基化的拟合关系，用于判断羟甲基化状态是否会影响TF的亲和力。
library(ggplot2)
library(ggrepel)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20210103.motif.analysis")
#1.全部TF做散点图
rt=read.csv("cor.mean.vaf.vs.score.allele.diff.csv",header = T)
rt$group="no.prefer"
sel1=which(rt$cor.p.value<0.05&rt$cor.estimate>0.7)
sel2=which(rt$cor.p.value<0.05&rt$cor.estimate< -0.7)

rt[sel1,]$group="prefer.alt.allele"
rt[sel2,]$group="prefer.ref.allele"
rt$logpvalue=-log10(rt$cor.p.value)
rt$absR2=abs(rt$cor.estimate)

p3=ggplot(rt,aes(y=absR2,x=logpvalue,color=group))+geom_vline(xintercept = -log10(0.05),color="red",linetype="dashed")+
    geom_point(size=3)+scale_color_manual(values = alpha(c("#666666","#3C5488","#E64B35"),0.9))+
    labs(y=expression(r^2),x="-log10 (Pvalue)",title="")+
    theme_classic(base_size = 15)+
    geom_text_repel(data=subset(rt,rt$cor.p.value<0.05),aes(label=TF),size=3, fontface="bold",force = T, segment.color = "black", show.legend = FALSE)

#2.只挑NRF1这个TF作例子 

cor_eqn=function(data){
    r=cor.test(data$mean.vaf,data$score.allele.diff,method = "pearson")
    eq <- substitute(italic(P.value) == p ~~ ","~~italic(r)^2~"="~rr,
                     list(p=as.numeric(format(r$p.value, digits = 3)),
                          rr=as.numeric(format(r$estimate, digits = 3))
                     ))
    as.expression(eq)
}					#粘贴在图中的公式

result=read.table("807.psy.ASH.motif.VAF.txt",head=T,sep="\t")
result$id2=paste(result$rsid,result$geneSymbol,sep="_")
result=result[!duplicated(result$id2),]
result$score.allele.diff=abs(result$scoreAlt-result$scoreRef)
tmp=result[result$geneSymbol=="NRF1",]

r=cor.test(tmp$mean.vaf,tmp$score.allele.diff,method = "pearson")
p=as.numeric(format(r$p.value, digits = 3))
rr=as.numeric(format(r$estimate, digits = 3))

p4=ggplot(tmp,aes(y=score.allele.diff,x=mean.vaf))+geom_point(size=3)+
  geom_smooth(method = "lm",color="black")+labs(y="binding likehood differences",x="VAF",title="")+
  theme_classic(base_size = 15)+annotate(x=0.2,y=9.9,label=cor_eqn(tmp),geom ="text" )
  
#3.饼图展示
dt=data.frame(table(rt$group))
dt[dt$Var1=="no.prefer",]$Freq=191+dt[dt$Var1=="no.prefer",]$Freq

p5=ggplot(dt, aes(x = "", y = Freq, fill = Var1)) +
  geom_bar(stat = "identity", width = 1)+
  coord_polar(theta = "y") + 
  labs(x = "", y = "", title = "") + 
  theme(axis.ticks = element_blank()) + 
  theme(legend.title = element_blank(), legend.position = "top") + 
  scale_fill_discrete(breaks = dt$Var1, labels = dt$Freq) + scale_fill_manual(values = alpha(c("#666666","#3C5488","#E64B35"),0.9))+theme_bw()+
  theme(
    panel.grid = element_blank(),
    panel.border= element_blank(),
    axis.text.y = element_blank(),
    axis.text.x = element_blank(),
    axis.ticks = element_blank(),
    axis.title = element_blank()
  )+guides(fill=F) 

layout <- matrix(c(1, 2, 3, 4,5,6), nrow = 3)
multiplot(plotlist=list(p1,p5,p4,p2,p3),layout=layout)